'''
Created on May 2, 2011

@author: oabalbin
'''
import sys
import numpy as np
from optparse import OptionParser
from collections import deque, defaultdict

def read_file(ifile, ofile):
    '''
    '''
    ifile = open(ifile)
    ofile =open(ofile,'w')
    
    for l in ifile:
        if l.startswith('#'):
            ofile.write(l)
            continue
        fields=l.strip('\n').replace(',','$').split('\t')
        
        pos= '%s@%s'%(fields[0],fields[1])
        #pos= '%s@%s|%s>%s'%(fields[0],fields[1],fields[3],fields[4])
        ofile.write(pos+',')
    
    ifile.close()
    ofile.close()
        

def make_posbase_file(ifile, ofile):
    '''
    Read a vcf file and write a pos.base file
    '''
    ifile = open(ifile)
    ofile =open(ofile.replace('.pos','.pos.base'),'w')
    
    for l in ifile:
        if l.startswith('#'):
            ofile.write(l)
            continue
        fields=l.strip('\n').replace(',','$').split('\t')
        
        #pos= '%s@%s'%(fields[0],fields[1])
        pos= '%s@%s|%s>%s'%(fields[0],fields[1],fields[3],fields[4])
        ofile.write(pos+',')
    
    ifile.close()
    ofile.close()

def read_position_base_file(ifile):
    '''
    '''
    ifile = open(ifile)
    for l in ifile:
        if l.startswith('#'):
            continue
        fields=l.strip('\n').split(',')
        fields.pop()
    
    ifile.close()
    
    positions, bases = deque(), deque()
    pos_dict=defaultdict()
    for snv in fields:
        #print snv
        pos,base=snv.split('|')[0],snv.split('|')[1]
        #positions.append(pos)
        #bases.append(base)
        pos_dict[pos]=base
    
    return pos_dict #positions, bases

    
            
def read_positions(ifile):
    '''
    '''
    ifile = open(ifile)
    for l in ifile:
        if l.startswith('#'):
            continue
        fields=l.strip('\n').split(',')
        fields.pop()
    
    ifile.close()
    return fields


def isec_arrays(filea,fileb):
    '''
    '''
    m1=set(read_positions(filea))
    m2=set(read_positions(fileb))
    #print sys.getsizeof(m1), sys.getsizeof(m2)
    isec=m1.intersection(m2)        
    return isec




def snv_concordance_to_database(query_snvs, ifile):
    '''
    common_snvs2databse: SNVs positions that are common with the database (set).
    snvs_in_database: ifile SNVs in the database in the format pos.base
    '''
    # read the pos.base file into two lists, the position list and the base list
    db_positions = read_position_base_file(ifile)
    # isec the commons_snvs2database and snvs_in_database
    isec_snvs = set(query_snvs.keys()).intersection(set(db_positions.keys()))
    # for the common snvs evaluate concordance.
    nrow = len(query_snvs)
    ncol = 1
    amat = np.zeros( (nrow,ncol) )
    
    for snv in isec_snvs:
        i = query_snvs[snv]
        db_base = db_positions[snv].split('>')[1]
        if db_base=='.':
            print db_base, snv
            amat[i,0] = 1
#    print amat.shape
#    print len(amat[amat==1])
    return np.reshape(amat,(nrow,-1))

if __name__ == '__main__':
    '''
    Write to an array the positions in a particualr vcf.
    
    '''
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-f", "--vcf_files", action="append", type="str",
                            dest="vcf_files",
                            help="others_vcf_files. Pass vcf files without .gz extention even if the file .gz is there")
    
    (options, args) = optionparser.parse_args()
    
    for ifile in options.vcf_files:
        ofile=ifile.replace('.vcf', '.pos')
        read_file(ifile, ofile)
        print >> sys.stderr,"DONE with file %s !"%(ifile)    
    
    
    
    
    '''
    for ifile in options.vcf_files:
        ofile=ifile.replace('.vcf', '.pos')
        make_posbase_file(ifile, ofile)
        print >> sys.stderr,"DONE with file %s !"%(ifile)
    
    query =  read_position_base_file(options.vcf_files[0])
    #snvs_in_database = read_position_base_file(options.vcf_files[1])
    query_snvs=defaultdict()
    for i, snv in enumerate(query):
        query_snvs[snv.split('|')[0]]=i 
        
    amat = snv_concordance_to_database(query_snvs, options.vcf_files[1])
    '''
    




